#### 参数设置 ####
projectPath <-  "D:\\work"     # 设置程序工作目录
cancer <-  "STAD"         # 设置访问的癌症名称
library(pacman)
##### Sept1.RNA表达数据下载与整理 #####
cancer.name <- paste("TCGA",cancer,sep = "-")  # 名称重定向
# 安装包
if('pacman' %in% .packages(all.available = TRUE)==FALSE){
  install.packages("pacman", repos = 'https://mirror.lzu.edu.cn/CRAN/')
}
library(pacman)
p_load(TCGAbiolinks, dplyr,tidyverse, magrittr, data.table, biomaRt)
# 设置工作目录
dataPath <-  paste(projectPath, "Data", sep = "\\")
if(!dir.exists(dataPath)) dir.create(dataPath)  # 新建一个Data文件夹，用于存放从TCGA上下载的数据
setwd(dataPath)
# 数据下载
for (data_type in c("HTSeq - Counts","HTSeq - FPKM-UQ")) {
  repeat{
    query.function <-  GDCquery(project = cancer.name,
                                legacy = FALSE, experimental.strategy = "RNA-Seq", 
                                data.category = "Transcriptome Profiling", data.type = "Gene Expression Quantification", 
                                workflow.type = data_type)
    if(class(query.function) != "try-error"){
      break
    }else{
      print("Try to connect to GDC...")
    }
  } # api访问函数
  GDCdownload(query.function, files.per.chunk = 100)   # 分批下载，断点续传
  # 数据预处理
  dataAssy = GDCprepare(query.function, summarizedExperiment = F)
  # 数据整理
  exp_data <-  dataAssy %>% 
    filter(str_detect(X1, "^ENSG")) %>% # 提取ENSG*注释的基因
    mutate(X1 = gsub("\\..*", "", X1)) %>% # 删除版本号
    dplyr::rename(ensembl_gene_id = X1) %>% # "X1"重命名为"ensembl_gene_id"
    set_colnames(str_replace(colnames(.), "(-\\w*){3}$", "")) # 修改TCGA样本名
  if(data_type == "HTSeq - FPKM-UQ"){
    # log2转化，如果是count值的话将直接输出，但是如果是FPKM值的话将要log2处理，使其服从正态分布
    exp_data[,2:ncol(exp_data)] = log2(exp_data[,2:ncol(exp_data)] + 1)
  }
  # 输出
  temp <- paste(cancer,"Portal_RNA_",sep = "_")
  outfile <-  paste0(temp, str_match(data_type, "- ([^\\s]*)$")[,2], ".txt")
  fwrite(exp_data, outfile, row.names = F, sep = "\t", quote = F)
}
##### Sept2.lncRNA注释 #####
p_load(clusterProfiler,biomaRt,dplyr,data.table,tibble,stringr,rtracklayer)
lncRNA_annotation <-  paste(projectPath, "lncRNA_annotation", sep = "\\")
if(!dir.exists(lncRNA_annotation)) dir.create(lncRNA_annotation)
temp <- paste(projectPath,"Data",sep = "\\")
temp <- paste(temp,cancer,sep = "\\")
temp <- paste(temp,"Portal_RNA_FPKM-UQ.txt",sep = "_")
setwd(projectPath)
exp_data_FPKM <-  fread(temp, h = T, sep = "\t", check.names = F) %>% column_to_rownames("ensembl_gene_id")
exp_data_T <-  exp_data_FPKM %>% dplyr::select(str_which(colnames(.), "-01A$")) # 提取癌症样本
nT <-  ncol(exp_data_T) # 373
exp_data_N <-  exp_data_FPKM %>% dplyr::select(ends_with("-11A"))
nN <-  ncol(exp_data_N) # 32
# 合并表达数据
exp_data_FPKM <-  cbind(exp_data_N, exp_data_T)
dim(exp_data_FPKM) # [1] 60483  405
exp_data <-  fread(temp, h = T, sep = "\t", check.names = F) %>% column_to_rownames("ensembl_gene_id")
exp_data_T <-  exp_data %>% dplyr::select(str_which(colnames(.), "-01A$")) # 提取癌症样本
nT <-  ncol(exp_data_T) # 373
exp_data_N <-  exp_data %>% dplyr::select(ends_with("-11A"))
nN <-  ncol(exp_data_N) # 32
# 合并表达数据
exp_data_final <-  cbind(exp_data_N, exp_data_T)
dim(exp_data_final) # [1] 60483  419
SampleLabel <-  factor(c(rep("NT", nN), rep(cancer, nT)), levels = c("NT", cancer))
rm(list = c("exp_data_N","exp_data_T"))
setwd(projectPath)
Anno_Data_lncRNA <-  import("gencode.v27.long_noncoding_RNAs.gtf.gz")
index <-  which(Anno_Data_lncRNA$type == 'gene')
length(index)
Tag <-  data.frame(Ensembl_ID = Anno_Data_lncRNA$gene_id[index],
                   Symbol = Anno_Data_lncRNA$gene_name[index],
                   Biotype = Anno_Data_lncRNA$gene_type[index])
Tag$Ensembl_ID <-  gsub('\\..*',"",Tag$Ensembl_ID)
Intersection_Set <-  intersect(Tag$Ensembl_ID,rownames(exp_data_final))
Intersection_Set <- as.data.frame(Intersection_Set)
exp_data_final <- exp_data_final[Intersection_Set$Intersection_Set,]
Tag <- Tag[match(rownames(exp_data_final),Tag$Ensembl_ID),]
tmp <- by(exp_data_final,
          Tag$Symbol,
          function(x) rownames(x)[which.max(rowMeans(x))])
prober <- as.character(tmp)
exp_data_final <- exp_data_final[rownames(exp_data_final) %in% prober,]
Tag <- Tag[match(rownames(exp_data_final),Tag$Ensembl_ID),]
rownames(exp_data_final) <- as.vector(as.matrix(Tag$Symbol))
setwd(lncRNA_annotation)
write.csv(exp_data_final,"FPKM_lncRNA_annotation.csv")
##### Sept3.RNA表达数据差异表达分析 #####
p_load(edgeR, ggbeeswarm, ggsignif,raster,rtracklayer)
deaPath <-  paste(projectPath, "DEA_NT", sep = "\\")
if(!dir.exists(deaPath)) dir.create(deaPath)  # 新建一个DEA_NT文件夹，用于存放差异表达分析结果
dataPath <-  paste(projectPath, "Data", sep = "\\")
temp <- paste(dataPath,cancer,sep = "\\")
temp <- paste(temp,"Portal_RNA_Counts.txt",sep = "_")
exp_data <-  fread(temp, h = T, sep = "\t", check.names = F) %>% column_to_rownames("ensembl_gene_id")
exp_data_T <-  exp_data %>% dplyr::select(str_which(colnames(.), "-01A$")) # 提取癌症样本
nT <-  ncol(exp_data_T) # 373
exp_data_N <-  exp_data %>% dplyr::select(ends_with("-11A"))
nN <-  ncol(exp_data_N) # 32
# 合并表达数据
exp_data_final <-  cbind(exp_data_N, exp_data_T)
dim(exp_data_final) # [1] 60483  419
SampleLabel <-  factor(c(rep("NT", nN), rep(cancer, nT)), levels = c("NT", cancer))
rm(list = c("exp_data_N","exp_data_T"))
setwd(projectPath)
Anno_Data_lncRNA <-  import("gencode.v27.long_noncoding_RNAs.gtf.gz")
index <-  which(Anno_Data_lncRNA$type == 'gene')
length(index)
Tag <-  data.frame(Ensembl_ID = Anno_Data_lncRNA$gene_id[index],
                   Symbol = Anno_Data_lncRNA$gene_name[index],
                   Biotype = Anno_Data_lncRNA$gene_type[index])
Tag$Ensembl_ID <-  gsub('\\..*',"",Tag$Ensembl_ID)
Intersection_Set <-  intersect(Tag$Ensembl_ID,rownames(exp_data_final))
Intersection_Set <- as.data.frame(Intersection_Set)
exp_data_final <- exp_data_final[Intersection_Set$Intersection_Set,]
Tag <- Tag[match(rownames(exp_data_final),Tag$Ensembl_ID),]
tmp <- by(exp_data_final,
          Tag$Symbol,
          function(x) rownames(x)[which.max(rowMeans(x))])
prober <- as.character(tmp)
exp_data_final <- exp_data_final[rownames(exp_data_final) %in% prober,]
Tag <- Tag[match(rownames(exp_data_final),Tag$Ensembl_ID),]
rownames(exp_data_final) <- as.vector(as.matrix(Tag$Symbol))
### edgeR差异表达分析 ###
# 构建DGEList对象
dge <-  DGEList(counts = exp_data_final, group = SampleLabel)
# filter
keep = rowSums(cpm(dge) > 0.5) > nN # 在至少较少组样本量中的表达量cpm值大于0.5
dge_filter = dge[keep, keep.lib.sizes = FALSE]
dim(dge_filter$counts) # [1] 13111   405
# TMM normalization：The observed counts should be scaled to the library sizes before comparing the read counts.
dge_norm = calcNormFactors(dge_filter)  # 默认 method = "TMM"
# 获取logCPM,用于可视化展示
logcpm = cpm(dge_norm, log = TRUE, prior.count = 1)
# Estimating the dispersion
design = model.matrix(~SampleLabel)
y = estimateDisp(dge_norm, design, robust = TRUE)
# DEA
fit = glmQLFit(y, design)
res = glmQLFTest(fit)
result = topTags(res, n = Inf)$table
result_sift <- result[result$FDR<0.01&abs(result$logFC)>2,]
# 输出
setwd(deaPath)
write_csv(data.frame(Symbol = rownames(result), result, check.names = F), "RNA_DEA_result_info.csv")
write_csv(data.frame(Symbol = rownames(result_sift), result_sift, check.names = F), "RNA_DEA_result_sift_info.csv")

##### Sept4.下载生存数据以及单因素cox分析与Lasso回归分析#####
Survive_Path <-  paste(projectPath, "Survive_Data", sep = "\\")
if(!dir.exists(Survive_Path)) dir.create(Survive_Path)  # 新建一个Survive_Data文件夹，用于存放从TCGA上下载的数据
setwd(Survive_Path)
COX_Path <- paste(Survive_Path,"COX",sep = "\\")
Lasso_Path <- paste(Survive_Path,"LASSO",sep = "\\")
dir.create(Lasso_Path)
dir.create(COX_Path)
p_load(data.table)
library(tidyverse)
library(UCSCXenaTools)
library(survival)
library(glmnet)
for (data_tpye in c("clinical","survival")) {
  describe_type <- XenaHub(hostName = "tcgaHub") %>%
    XenaFilter(filterDatasets = data_tpye) %>%
    XenaFilter(filterDatasets = cancer)
  data_query <- XenaQuery(describe_type)
  data_dowload <- XenaDownload(data_query)
  if (data_tpye == "clinical") {
    cilinical <- XenaPrepare(data_dowload)
    write.csv(cilinical,"cilinical.csv")
  }else{
    survival <- XenaPrepare(data_dowload)%>%
      filter(OS.time >0 & !is.na(OS) ) %>%
      droplevels
    survival <- survival[,c(2,3,4)]
    write.csv(survival,"survival.csv")
  }
}
deaPath <-  paste(projectPath, "DEA_NT", sep = "\\")
lncRNA_annotation <-  paste(projectPath, "lncRNA_annotation", sep = "\\")
setwd(deaPath)
gene_ex_sift <- read.csv("RNA_DEA_result_sift_info.csv",fill = T,sep = ",")
setwd(lncRNA_annotation)
gene <- read.csv("FPKM_lncRNA_annotation.csv",fill = T,sep = ",")
gene <- gene[gene$X %in% gene_ex_sift$Symbol,]
gene <- t(gene)
colnames(gene) <- gene[1,]
gene <- gene[-1,]
rownames(gene) <- gsub("\\.","\\-",rownames(gene))
gene <- gene[grepl("*-01A",rownames(gene)),]
rownames(gene) <- substr(rownames(gene),1,12)
write.table(gene,file = "lnc_gene.txt",sep = "\t",quote=F)
gene <- fread("lnc_gene.txt", sep = "\t")
ex_cl <- merge(survival,gene,by.x = "_PATIENT",by.y = "V1")
ex_cl <- unique(ex_cl)
rownames(ex_cl) <- ex_cl[,1]
ex_cl <- ex_cl[,-1]
coxR <- data.frame()
coxf <- function(x){
  fmla1 <- as.formula(Surv(OS.time,OS)~ex_cl[,x])
  mycox <- coxph(fmla1,data=ex_cl)
}
for(a in colnames(ex_cl[,3:ncol(ex_cl)])){
  mycox <- coxf(a)
  coxResult <-  summary(mycox)
  coxR <- rbind(coxR,cbind(lncname = a,
                           HR = coxResult$coefficients[,"exp(coef)"],
                           P = coxResult$coefficients[,"Pr(>|z|)"]))
}
setwd(COX_Path)
coxR <- coxR[coxR$P < 0.05,]
write.table(coxR,"coxResult.txt",sep = "\t",quote = F,row.names=F)
ex_cl$OS.time <- ex_cl$OS.time/365
v1 <- as.matrix(ex_cl[,c(3:ncol(ex_cl))])
v2 <- as.matrix(Surv(ex_cl$OS.time,ex_cl$OS))
myfit <- glmnet(v1, v2, family = "cox")
setwd(Lasso_Path)
pdf("lambda.pdf")
plot(myfit, xvar = "lambda", label = TRUE)
dev.off()
myfit2 <- cv.glmnet(v1, v2, family="cox")
pdf("min.pdf")
plot(myfit2)
abline(v=log(c(myfit2$lambda.min,myfit2$lambda.1se)),lty="dashed")
dev.off()
myfit2$lambda.min
coe <- coef(myfit, s = myfit2$lambda.min)
act_index <- which(coe != 0)
act_coe <- coe[act_index]
result <- row.names(coe)[act_index]
write.table(result,"Lasso_result.txt",quote = F,sep = "\t")


##### Sept5.多因素cox回归分析#####
Multifactor_Path <-  paste(Survive_Path, "Multifactor", sep = "\\")
setwd(Survive_Path)
dir.create(Multifactor_Path)
setwd(COX_Path)
Single_factor <- read.table("coxResult.txt",sep = "\t",header = T)
setwd(Lasso_Path)
Lasso <- read.table("Lasso_result.txt",sep = "\t",header = T)
index <- merge(Lasso,Single_factor,by.x = "x",by.y = "lncname")
setwd(Survive_Path)
survival <- read.csv("survival.csv",sep = ",",fill = T,check.names = T)
survival <- survival[,-1]
setwd(lncRNA_annotation)
gene <- read.csv("FPKM_lncRNA_annotation.csv",fill = T,sep = ",")
gene <- gene[gene$X %in% index$x,]
gene <- t(gene)
colnames(gene) <- gene[1,]
gene <- gene[-1,]
rownames(gene) <- gsub("\\.","\\-",rownames(gene))
gene <- gene[grepl("*-01A",rownames(gene)),]
rownames(gene) <- substr(rownames(gene),1,12)
write.table(gene,file = "lnc_gene.txt",sep = "\t",quote=F)
gene <- fread("lnc_coxre.txt", sep = "\t")
file.remove("lnc_coxre.txt")
ex_cl <- merge(survival,gene,by.x = "X_PATIENT",by.y = "V1")
ex_cl <- unique(ex_cl)
rownames(ex_cl) <- ex_cl[,1]
ex_cl <- ex_cl[,-1]
ex_cl$OS.time <- ex_cl$OS.time/365
fmla1 <- as.formula(Surv(OS.time,OS)~.)
mycox <- coxph(fmla1,data = ex_cl)
risk_score <- predict(mycox,type = "risk",newdata = ex_cl)
risk_level <- as.factor(ifelse(risk_score > median(risk_score),"High","Low"))
setwd(Multifactor_Path)
write.table(cbind(id=rownames(cbind(ex_cl[,1:2],
                                    risk_score,risk_level)),
                  cbind(ex_cl[,1:2],risk_score,risk_level)),
            "risk_score.txt",sep="\t",quote=F,row.names=F)
summary(mycox)
library(survminer)
pdf("forest1.pdf",12,8)
ggforest(mycox,fontsize = 1)
dev.off()
library(ggplot2)
data <- read.table("risk_score.txt",
                   header = T,sep = "\t",check.names = F)
attach(data)
pdf("risk socre in scatter diagram.pdf")
plot(risk_score,OS.time,
     col = ifelse(OS == 1,"#F0EC90","#77D980"),
     xlab = "Patients (increasing risk socre)",ylab="Survival time(years)",pch=16)
legend("topright", c("Death", "Alive"), pch=16, col=c("#F0EC90","#77D980"))
dev.off()
setwd(deaPath)
de_set <- read.csv("RNA_DEA_result_info.csv",fill = T,sep = ",")
df <- de_set
rownames(df) <- df[,1]
df <- df[,-1]
df$threshold = factor(ifelse(df$FDR < 0.05 & abs(df$logFC) >= 2, 
                             ifelse(df$logFC>= 2 ,'Up','Down'),'NoSignifi'),
                      levels=c('Up','Down','NoSignifi'))
deg <- df[rownames(df) %in% index$x,]

p <- ggplot(df,aes(x=logFC,y= -log2(PValue),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("#F0EC90","#77D980","#90E5F0"))+#确定点的颜色
  theme_bw()+#修改图片背景
  theme(
    legend.title = element_blank()#不显示图例标题
  )+
  ylab('-log10 (p-Value)')+#修改y轴名称
  xlab('log(FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-1,1),lty=3,col="black",lwd=0.5) +
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=0.5)

k<- p + geom_point(size = 3, shape = 1, data = deg) +
  ggrepel::geom_label_repel(
    aes(label = rownames(deg)),
    size=5,
    data = deg,
    color="#7E77D9"
  )
pdf("Volcano.pdf")
k
dev.off()





##### Sept6.ROC#####
p_load(timeROC,survival)
ROC_Path <-  paste(projectPath, "ROC", sep = "\\")
setwd(projectPath)
dir.create(ROC_Path)
setwd(Multifactor_Path)
risk_score <- read.table("risk_score.txt",header=T,sep="\t")
predict_3_year <- 3
predict_5_year <- 5
ROC<-timeROC(T = risk_score$OS.time,delta = risk_score$OS,
             marker = risk_score$risk_score , cause=1,
             weighting = "marginal",
             times = c(predict_3_year,predict_5_year),ROC=TRUE)
setwd(ROC_Path)
pdf("ROC.pdf")
plot(ROC,time = predict_3_year,col="#F0EC90",title=FALSE,lwd = 3)
plot(ROC,time = predict_5_year,col="#51F588",add = TRUE,title = FALSE,lwd = 3)
legend("bottomright",
       c(paste("AUC of 3 year survival: ",round(ROC$AUC[1],3)),
         paste("AUC of 5 year survival: ",round(ROC$AUC[2],3))),col = c("#F0EC90","#51F588"),lwd = 3)
dev.off()



##### Sept6.浸润分析 #####
p_load(data.table,stringr)
CibersortPath <-  paste(projectPath, "Cibersort_Analyse", sep = "\\")
if(!dir.exists(CibersortPath)) dir.create(CibersortPath)
temp <- paste(projectPath,"Data",sep = "\\")
temp <- paste(temp,cancer,sep = "\\")
temp <- paste(temp,"Portal_RNA_FPKM-UQ.txt",sep = "_")
setwd(projectPath)
exp_data_FPKM <-  fread(temp, h = T, sep = "\t", check.names = F) %>% tibble::column_to_rownames("ensembl_gene_id")
exp_data_T <-  exp_data_FPKM %>% dplyr::select(str_which(colnames(.), "-01A$")) # 提取癌症样本
nT <-  ncol(exp_data_T) # 373
exp_data_N <-  exp_data_FPKM %>% dplyr::select(ends_with("-11A"))
nN <-  ncol(exp_data_N) # 32
# 合并表达数据
exp_data_FPKM <-  cbind(exp_data_N, exp_data_T)
dim(exp_data_FPKM) # [1] 60483  405
source("Cibersort.R")
p_load(clusterProfiler,biomaRt)
# biomaRt包访问ensembl数据库
ensembl <-  useEnsembl("ensembl", dataset = "hsapiens_gene_ensembl")
genes_info <-  getBM(attributes = c("ensembl_gene_id", "external_gene_name"), filters = "ensembl_gene_id", values = rownames(exp_data_FPKM), mart = ensembl) # 无对应Symbol的自动删除
genes_info <- genes_info[-which(genes_info$external_gene_name==""),]
exp_data_FPKM <- exp_data_FPKM[rownames(exp_data_FPKM) %in% genes_info$ensembl_gene_id,]
genes_info <- genes_info[match(rownames(exp_data_FPKM),genes_info$ensembl_gene_id),]
tmp <- by(exp_data_FPKM,
          genes_info$external_gene_name,
          function(x) rownames(x)[which.max(rowMeans(x))])
prober <- as.character(tmp)
exp_data_FPKM <- exp_data_FPKM[rownames(exp_data_FPKM) %in% prober,]
genes_info <- genes_info[match(rownames(exp_data_FPKM),genes_info$ensembl_gene_id),]
exp_data_FPKM_L <- exp_data_FPKM
rownames(exp_data_FPKM) <- genes_info$external_gene_name
exp_data_FPKM <- rbind(ID=colnames(exp_data_FPKM),exp_data_FPKM)
write.table(exp_data_FPKM,file="ready.txt",sep="\t",quote=F,col.names=F)
results=CIBERSORT("LM22.txt", "ready.txt", perm=1000, QN=TRUE)
file.remove("ready.txt")
setwd(CibersortPath)
results <- results[results[,23] < 0.05,] #233
results <- results[,c(1:22)]
setwd(CibersortPath)
write.csv(results,"Cibersort_results.csv")
############绘图区
p_load(RColorBrewer,dplyr,tidyr,ggplot2)
library(plyr)
####正常人与患者的免疫细胞水平
if(T){
  mypalette <- colorRampPalette(brewer.pal(8,"Set3"))
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat$Sample <- substr(x = dat$Sample,start = 14,stop = length(dat$Sample))
  dat$Sample <- gsub("*01A","tumor",dat$Sample)
  dat$Sample <- gsub("*11A","normal",dat$Sample)
  ggplot(dat,aes(x = Sample,y = Proportion,fill = Cell_type)) +
    geom_bar(stat = "identity") +
    labs(fill = "Cell Type",y = "Estiamted Proportion",title = "The normal and tumor person of immune cells in 22") +
    theme_bw() +
    scale_y_continuous(expand = c(0.01,0)) +
    scale_fill_manual(values = mypalette(22))+
    theme(plot.title = element_text(hjust = 0.5))
  ggsave('The difference between normal and abnormal of immune cells in 22.png',width = 400,height = 230
         ,units = "mm")
  
}
#####正常人与患者的免疫细胞水平,以百分比形式
if(T){
  p_load(plyr)
  mypalette <- colorRampPalette(brewer.pal(8,"Set3"))
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat$Sample <- substr(x = dat$Sample,start = 14,stop = length(dat$Sample))
  dat$Sample <- gsub("*01A","tumor",dat$Sample)
  dat$Sample <- gsub("*11A","normal",dat$Sample)
  date <- ddply(dat, "Sample", transform, percent=Proportion / sum(Proportion))
  ggplot(date,aes(x = Sample,y = 100*percent,fill = Cell_type)) +
    geom_bar(stat = "identity") +
    labs(fill = "Cell Type",y = "%",title = "The normal and tumor person of immune cells in 22") +
    theme_bw() +
    scale_y_continuous(expand = c(0.01,0)) +
    scale_fill_manual(values = mypalette(22))+
    theme(plot.title = element_text(hjust = 0.5))
  ggsave('The difference percent between normal and abnormal of immune cells in 22.png',width = 400,height = 230
         ,units = "mm")
}
#####患者的免疫细胞占比
if(T){
  mypalette <- colorRampPalette(brewer.pal(8,"Set3"))
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat <- dat[grep("*01A",dat$Sample),]
  date <- ddply(dat, "Sample", transform, percent=Proportion / sum(Proportion))
  ggplot(date,aes(x = Sample,y = 100*percent,fill = Cell_type)) +
    geom_bar(stat = "identity") +
    labs(fill = "Cell Type",y = "%",title = "Tumor tissues") +
    theme_bw() +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          legend.position = "bottom") + 
    scale_y_continuous(expand = c(0.01,0)) +
    scale_fill_manual(values = mypalette(22))+
    theme(plot.title = element_text(hjust = 0.5))
  ggsave('Percentage of immune cells in patients.png',width = 400,height = 230
         ,units = "mm")
  
}
#####正常人的免疫细胞占比
if(T){
  mypalette <- colorRampPalette(brewer.pal(8,"Set3"))
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat <- dat[grep("*11A",dat$Sample),]
  date <- ddply(dat, "Sample", transform, percent=Proportion / sum(Proportion))
  ggplot(date,aes(x = Sample,y = 100*percent,fill = Cell_type)) +
    geom_bar(stat = "identity") +
    labs(fill = "Cell Type",y = "%",title = "Normal tissues") +
    theme_bw() +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          legend.position = "bottom") + 
    scale_y_continuous(expand = c(0.01,0)) +
    scale_fill_manual(values = mypalette(22))+
    theme(plot.title = element_text(hjust = 0.5))
  ggsave('Percentage of immune cells in normal.png',width = 400,height = 230
         ,units = "mm")
  
}
# ####小提琴图未完成
# if(T){
#   dat <- results %>%
#     as.data.frame() %>%
#     tibble::rownames_to_column("Sample") %>%
#     gather(key = Cell_type,value = Proportion,-Sample)
#   dat$Sample <- substr(x = dat$Sample,start = 14,stop = length(dat$Sample))
#   dat$Sample <- gsub("*01A","tumor",dat$Sample)
#   dat$Sample <- gsub("*11A","normal",dat$Sample)
#   # dat$Proportion <- dat$Proportion*100
#   ggviolin(dat, x = "Cell_type", y = "Proportion", color = "Sample")
#   
#   # 
#   # ggplot(dat,aes(x = Sample,y = Proportion,fill = Cell_type))+
#   #   geom_violin(alpha = 1,              # 透明度
#   #               trim = T,               # 是否修剪尾巴，即将数据控制到真实的数据范围内
#   #               scale = "width")        # 如果“area”(默认)，所有小提琴都有相同的面积(在修剪尾巴之前)。如果是“count”，区域与观测的数量成比例。如果是“width”，所有的小提琴都有相同的最大宽度。
#   ggsave('volin.png')
#   }
####总T细胞
if(T){
  mypalette <- colorRampPalette(brewer.pal(8,"Set3"))
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat$Sample <- substr(x = dat$Sample,start = 14,stop = length(dat$Sample))
  dat$Sample <- gsub("*01A","tumor",dat$Sample)
  dat$Sample <- gsub("*11A","normal",dat$Sample)
  dat <- dat[grep("T.cell*",dat$Cell_type),]
  ggplot(dat,aes(x = Sample,y = Proportion,fill = Cell_type)) +
    geom_bar(stat = "identity") +
    labs(fill = "Cell Type",y = "Immunocyte content",title = "All T cells") +
    theme_bw() +
    scale_y_continuous(expand = c(0.01,0)) +
    scale_fill_manual(values = mypalette(22))+
    theme(plot.title = element_text(hjust = 0.5))
  ggsave('All T cells.png',width = 400,height = 230
         ,units = "mm")
}
####总B细胞
if(T){
  mypalette <- colorRampPalette(brewer.pal(8,"Set3"))
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat$Sample <- substr(x = dat$Sample,start = 14,stop = length(dat$Sample))
  dat$Sample <- gsub("*01A","tumor",dat$Sample)
  dat$Sample <- gsub("*11A","normal",dat$Sample)
  dat <- dat[grep("B.cells.*",dat$Cell_type),]
  ggplot(dat,aes(x = Sample,y = Proportion,fill = Cell_type)) +
    geom_bar(stat = "identity") +
    labs(fill = "Cell Type",y = "Immunocyte content",title = "All B cells") +
    theme_bw() +
    scale_y_continuous(expand = c(0.01,0)) +
    scale_fill_manual(values = mypalette(22))+
    theme(plot.title = element_text(hjust = 0.5))
  ggsave('All B cells.png',width = 400,height = 230
         ,units = "mm")
}
####总巨噬细胞
if(T){
  mypalette <- colorRampPalette(brewer.pal(8,"Set3"))
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat$Sample <- substr(x = dat$Sample,start = 14,stop = length(dat$Sample))
  dat$Sample <- gsub("*01A","tumor",dat$Sample)
  dat$Sample <- gsub("*11A","normal",dat$Sample)
  dat <- dat[grep("Macrophages.*",dat$Cell_type),]
  ggplot(dat,aes(x = Sample,y = Proportion,fill = Cell_type)) +
    geom_bar(stat = "identity") +
    labs(fill = "Cell Type",y = "Immunocyte content",title = "All Macrophages cells") +
    theme_bw() +
    scale_y_continuous(expand = c(0.01,0)) +
    scale_fill_manual(values = mypalette(22))+
    theme(plot.title = element_text(hjust = 0.5))
  ggsave('All Macrophages cells.png',width = 400,height = 230
         ,units = "mm")
}
####箱式图
if(T){
  p_load(ggpubr)
  dat <- results %>%
    as.data.frame() %>%
    tibble::rownames_to_column("Sample") %>%
    gather(key = Cell_type,value = Proportion,-Sample)
  dat$Sample <- substr(x = dat$Sample,start = 14,stop = length(dat$Sample))
  dat$Sample <- gsub("*01A","tumor",dat$Sample)
  dat$Sample <- gsub("*11A","normal",dat$Sample)
  ggplot(dat,aes(x = Cell_type,y = Proportion,fill = Sample),
         add = "jitter",facet.by = "Sample")+
    geom_boxplot()+
    theme_bw()+
    stat_compare_means(aes(group = Sample,label=paste0("p=",..p.format..)),
                       size = 2.5,
                       method = "t.test",label.y = 0.3)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
  ggsave('All Immunize Cells.png',width = 400,height = 230
         ,units = "mm")
}
##### Sept7.相关性分析 #####
cor_Path <- paste(projectPath, "COR", sep = "\\")
setwd(projectPath)
dir.create(cor_Path)
setwd(COX_Path)
Single_factor <- read.table("coxResult.txt",sep = "\t",header = T)
setwd(Lasso_Path)
Lasso <- read.table("Lasso_result.txt",sep = "\t",header = T)
index <- merge(Lasso,Single_factor,by.x = "x",by.y = "lncname")
setwd(Survive_Path)
survival <- read.csv("survival.csv",sep = ",",fill = T,check.names = T)
survival <- survival[,-1]
setwd(lncRNA_annotation)
gene <- read.csv("FPKM_lncRNA_annotation.csv",fill = T,sep = ",")
gene <- gene[gene$X %in% index$x,]
gene <- t(gene)
colnames(gene) <- gene[1,]
gene <- gene[-1,]
rownames(gene) <- gsub("\\.","\\-",rownames(gene))
setwd(CibersortPath)
cibersoft <- read.csv("Cibersort_results.csv",fill = T,sep = ",")
index <- cibersoft$X
gene <- gene[index,]
write.csv(gene,"gene.csv")
gene <- read.csv("gene.csv",fill = T,sep = ",")
cor_ex <- merge(cibersoft,gene,by = "X")
rownames(cor_ex) <- cor_ex[,1]
cor_ex <- cor_ex[,-1]
cor_matrix <- cor(cor_ex,method = "spearman")
setwd(cor_Path)
library(corrplot) 
pdf("cor.pdf")
corrplot(cor_matrix,method = "color")
dev.off()
